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ABSTRACT 

The  direct  simulation  Monte  Carlo  method  (DSMC)  has  evolved  over  40  years  into  a  powerful  numerical 
technique  for  the  computation  of  complex,  nonequilibrium  gas  flows.  In  this  context,  nonequilibrium  means 
that  the  velocity  distribution  function  is  not  in  an  equilibrium  form  due  to  a  low  number  of  intermolecular 
collisions  within  a  fluid  element.  In  atmospheric  entry,  nonequilibrium  conditions  occur  at  high  altitude  and 
in  regions  of  flow  fields  with  small  length  scales.  In  this  first  article  of  two  parts,  the  theoretical  basis  of  the 
DSMC  technique  is  discussed.  In  addition,  the  methods  used  in  DSMC  are  described  for  simulation  of  high 
temperature,  real  gas  effects  and  gas-surface  interactions.  The  current  status  of  the  various  models  is 
reviewed  and  areas  where  further  work  is  required  are  identified. 


1.0  INTRODUCTION 


The  analysis  of  dilute  gas  flows  at  all  Knudsen  numbers  can  be  derived  from  the  Boltzmann  equation  that 
describes  the  evolution  of  the  molecular  velocity  distribution  function  (VDF)  [1],  In  the  absence  of  a  body 
force,  the  Boltzmann  equation  is  written: 

^(nf)+C.^(nf)  =  A(f)  (1.1) 

where  /  is  the  VDF,  n  is  the  number  density,  C  is  the  particle  velocity  vector,  r  is  the  particle  position 
vector,  t  is  time,  and  A  (A  represents  the  rate  of  change  in  the  VDF  due  to  collision  processes.  The 
equilibrium  solution  of  the  Boltzmann  equation  is  the  Maxwellian  VDF: 
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where  m  is  the  mass  of  a  particle,  k  is  Boltzmann’s  constant,  and  T  is  the  temperature.  The  physical 
mechanism  that  maintains  the  VDF  in  equilibrium  is  inter-molecular  collisions,  and  so  a  gas  falls  into  a 
nonequilibrium  state  under  conditions  where  there  is  not  a  large  enough  number  of  collisions  occurring  to 
maintain  equilibrium.  The  two  main  physical  flow  conditions  that  lead  to  nonequilibrium  are  low  density  and 
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small  length  scales.  A  low  density  leads  to  a  reduced  collision  rate  while  a  small  length  scale  reduces  the  size 
of  a  fluid  element.  The  usual  metric  for  determining  whether  a  particular  gas  flow  is  in  a  state  of 
nonequilibrium  is  the  Knudsen  number  defined  as  follows: 

Kn  =  —  (1.3) 

L 

where  A  is  the  mean  free  path  of  the  gas  and  L  is  the  characteristic  length  scale.  The  mean  free  path  is  the 
average  distance  traveled  by  each  particle  between  collisions  and  is  given  for  a  hard  sphere  by 
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where  n  is  again  the  number  density,  and  a  is  the  hard  sphere  collision  cross  section.  Thus,  at  low  density,  the 
mean  free  path  (and  therefore  Kn)  becomes  large.  Similarly,  for  small  length  scales,  L  becomes  small  and 
again  Kn  becomes  large.  As  a  guiding  rule,  it  is  generally  accepted  that  kinetic  nonequilibrium  effects 
become  important  when  Kn  >  0.01 . 


At  a  Knudsen  number  near  zero,  the  velocity  distribution  function  everywhere  in  a  flow  field  has  the 
Maxwellian  form,  there  are  no  molecular  transport  processes,  such  as  viscosity  and  thermal  conductivity,  and 
the  flow  may  be  modeled  using  the  Euler  equations.  Indeed,  the  Euler  equations  of  fluid  flow  can  be  derived 
by  taking  moments  of  the  Boltzmann  equation  and  evaluating  them  using  the  Maxwellian  VDF.  As  the 
Knudsen  number  increases  up  to  values  below  0.01,  the  velocity  distribution  function  in  the  flow  field  may  be 
represented  as  a  small  perturbation  from  the  equilibrium  Maxwellian  form  that  is  known  as  the  Chapman- 
Enskog  distribution  [1],  Evaluation  of  moments  of  the  Boltzmann  equation  using  the  Chapman-Enskog  VDF 
leads  to  the  Navier-Stokes  equations  in  which  shear  stress  and  heat  flux  depend  linearly  on  the  spatial 
gradients  of  velocity  and  temperature,  respectively.  As  the  Knudsen  number  rises  above  0.01,  these  linear 
transport  relations  are  unable  to  accurately  describe  the  strong  nonequilibrium  processes.  It  is  then  necessary 
to  develop  higher  order  sets  of  partial  differential  equations  (such  as  the  Burnett  equations)  or  to  solve  the 
Boltzmann  equation.  While  there  has  been  some  success  achieved  in  formulating  and  solving  the  Burnett 
equations,  there  remain  issues  with  boundary  conditions,  and  it  is  not  clear  that  the  amount  of  additional 
Knudsen  number  range  provided  is  worth  the  significant  additional  effort  in  numerical  analysis. 
Unfortunately,  development  of  robust  and  general  numerical  solution  schemes  for  the  Boltzmann  equation  has 
also  proved  a  significant  challenge.  Again,  some  progress  has  been  made,  but  there  is  still  much  work  to  be 
done  to  be  able  to  simulate  all  of  the  flow  physics  of  interest  in  real  world  problems. 


The  direct  simulation  Monte  Carlo  (DSMC)  method  was  first  introduced  by  Graeme  Bird  in  1961  [2] 
specifically  to  analyze  high  Knudsen  number  flows.  Since  that  time,  Bird  has  written  two  books  on  the 
method  [3,4]  and  thousands  of  research  papers  have  been  published  that  report  on  development  and 
application  of  the  technique.  The  significance  of  the  DSMC  technique  has  been  its  ability  over  40  years  of 
development  to  fill  the  void  described  above  in  gas  analysis  methodology  for  high  Knudsen  number  flows. 
The  DSMC  technique  emulates  the  same  physics  as  the  Boltzmann  equation  without  providing  a  direct 
solution.  The  DSMC  method  follows  a  representative  set  of  particles  as  they  collide  and  move  in  physical 
space.  It  has  been  demonstrated  that  DSMC  converges  to  solution  of  the  Boltzmann  equation  in  the  limit  of  a 
very  large  number  of  particles  [4], 

High  Knudsen  number  conditions  arise  in  many  areas  of  science  and  technology  including  space  and 
atmospheric  science,  vapor  processing  of  materials,  spacecraft  propulsion  systems,  and  micro-scale  gas  flows. 
Atmospheric  entry  flow  conditions  may  fall  into  the  kinetic  nonequilibrium  regime  at  sufficiently  low  density 
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(that  occurs  at  high  altitude  in  a  planet’s  atmosphere)  and  for  very  small  entry  shapes  (e.g.  meteoroids  that 
have  a  diameter  on  the  order  of  a  centimeter  [5]).  In  addition,  situations  arise  where  localized  regions  of  a 
flow  may  contain  low  density  (e.g.  the  wake  behind  a  capsule)  or  small  length  scales  (e.g.  sharp  leading  edges 
on  a  vehicle,  or  shock  waves  and  boundary  layers  that  may  have  very  steep  spatial  gradients  in  flow  field 
properties). 

For  hypersonic  flows,  it  is  a  valid  question  to  ask  whether  high  Knudsen  number  phenomena  lead  to  effects 
that  are  of  practical  significance.  Recent  detailed  studies  [6,7]  have  compared  DSMC  and  CFD  computations 
of  hypersonic  flows  over  cylinders  and  wedges  for  global  Knudsen  numbers  ranging  from  continuum 
(Kn=0.002)  to  rarefied  (Kn=0.25).  The  focus  of  the  studies  was  the  effect  of  any  nonequilibrium  flow 
phenomena  on  surface  properties  such  as  drag  and  heat  transfer.  It  was  found  at  Kn=0.002  that  DSMC  and 
CFD  gave  identical  results  for  all  surface  properties  including  drag  force  and  peak  heat  transfer.  However,  as 
Kn  was  increased,  the  differences  between  DSMC  and  CFD  surface  results  also  grew  larger.  For  example,  in 
Mach  25  flow  of  nitrogen  over  a  cylinder  at  the  highest  Knudsen  number  of  0.25,  in  comparison  to  CFD, 
DSMC  predicted  a  23%  lower  drag  force  and  a  29%  lower  peak  heat  flux.  These  differences  are  clearly 
significant  and  indicate  that  accurate  determination  of  surface  properties  in  high  Knudsen  number  flows  does 
require  a  non-continuum,  kinetic  approach,  such  as  DSMC. 

An  important  part  of  the  success  of  the  DSMC  technique  in  analyzing  high  Knudsen  number  atmospheric 
entry  flows,  has  been  the  ability  to  include  in  the  technique  models  that  are  effective  in  simulating  high 
temperature,  real  gas  effects.  Such  effects  include  mixtures  of  chemical  species,  relaxation  of  internal  energy 
modes,  chemical  reactions  such  as  dissociation  and  ionization,  radiation,  and  gas-surface  interaction.  In  this 
article,  the  fundamental  aspects  of  the  DSMC  technique  are  described  with  an  emphasis  on  physical  modeling 
issues  related  to  its  application  to  hypersonic,  atmospheric  entry  problems  and  the  simulation  of  the  associated 
real  gas  effects.  Examples  are  provided  that  illustrate  the  current  capabilities  of  these  models  and  areas  where 
further  work  is  needed  are  identified.  The  second  part  of  this  article  [8]  reviews  the  current  status  of  existing 
DSMC  codes  and  their  application  to  analysis  of  hypersonic  entry  flows. 


2.0  BASIC  ALGORITHM  OF  THE  DSMC  TECHNIQUE 

The  DSMC  technique  emulates  the  physics  of  the  Boltzmann  equation  by  following  the  motions  and  collisions 
of  a  large  number  of  model  particles.  Each  particle  possesses  molecular  level  information  including  a  position 
vector,  a  velocity  vector,  and  physical  information  such  as  mass  and  size.  Particle  motion  and  collisions  are 
decoupled  over  a  time  step  At  that  is  smaller  than  the  local  mean  free  time.  During  the  movement  of  particles, 
boundary  conditions  such  as  reflection  from  solid  surfaces  are  applied.  The  physical  domain  to  be  simulated 
in  a  DSMC  computation  is  covered  by  a  mesh  of  cells.  These  cells  are  used  to  collect  together  particles  that 
may  collide.  There  are  a  number  of  DSMC  schemes  for  simulating  collisions  and  all  of  them  achieve  a  faster 
numerical  performance  than  the  molecular  dynamics  (MD)  method  [9]  by  ignoring  the  influence  of  the 
relative  positions  of  particles  within  a  cell  in  determining  particles  that  collide.  This  simplification  requires 
that  the  size  of  each  cell  be  less  than  the  local  mean  free  path  of  the  flow.  Bird’s  No  Time  Counter  (NTC) 
scheme  [4]  is  the  most  widely  used  collision  scheme  in  which  a  number  of  particle  pairs  in  a  cell  are  formed 
that  is  given  by: 


N=  —  nN(og)  At 

c  ^  n  ‘-’/max 


(2.1) 
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where  n  is  the  number  density,  N  is  the  average  number  of  particles  in  the  cell,  cr  is  the  collision  cross 
section,  and  g  is  the  relative  velocity.  Each  of  the  Nc  pairs  of  particles  is  formed  at  random  regardless  of 
position  in  the  cell,  and  then  a  probability  of  collision  for  each  pair  is  evaluated  using: 


P  = 


(2.2) 


This  procedure  reproduces  the  expected  equilibrium  collision  rate  under  conditions  of  equilibrium.  It  is 
determined  whether  the  particle  pair  actually  collides  by  comparing  the  collision  probability,  Pc,  to  a  random 
number.  When  a  collision  occurs,  post-collision  velocities  are  calculated  using  conservation  of  momentum 
and  energy. 


The  cells  employed  for  simulating  collisions  are  also  often  used  for  the  sampling  of  macroscopic  flow 
properties  such  as  density,  velocity,  and  temperature.  There  is  no  necessity  to  have  the  collision  and  sampling 
cells  be  identical,  however,  and  sometimes  a  coarser  mesh  is  used  for  sampling. 


The  basic  steps  in  each  iteration  of  the  DSMC  method  are:  (1)  move  particles  over  the  time  step  At;  (2)  apply 
boundary  conditions  such  as  introducing  new  particles  at  inflow  boundaries,  removing  particles  at  outflow 
boundaries,  and  processing  reflections  at  solid  boundaries;  (3)  sort  particles  into  cells  and  calculate  collisions; 
and  (4)  sample  average  particle  information.  As  an  example  of  how  sampled  particle  information  is  employed 
to  determine  a  macroscopic  flow  property,  the  average  mass  density  in  a  computational  cell  of  volume  V  is 
given  by 


P  = 


i=i 

VxN, 


(2.3) 


where  m,  is  the  mass  of  particle  i,  Np  is  the  total  number  of  particles  that  have  occupied  this  cell  over  N, 
iterations  of  the  computation.  Also,  as  an  illustration  of  the  determination  of  surface  properties,  for  a  surface 
element  of  area  A,  the  average  shear  stress  is  given  by 


-  ut)j 

j _ 

AxNtxAt 


(2.4) 


where  u\  and  u\  are  the  incident  and  reflected  velocity  components  tangent  to  the  surface  element  of  each 
particle  j  to  hit  the  element  during  N,  iterations  of  the  computation. 


A  simulation  begins  from  some  initial  condition,  and  a  finite  number  of  iterations  must  elapse  in  order  for  the 
flow  to  reach  a  steady  state.  Generally,  steady  state  is  detected  as  a  leveling  off  of  the  total  number  of 
particles  in  the  simulation.  After  steady  state  is  reached,  sampling  of  flow  field  and  surface  properties  begins 
and  the  simulation  is  continued  a  further  number  of  iterations  in  order  to  reduce  the  statistical  noise  in  the 
sampled  information  to  an  acceptable  level.  A  typical  DSMC  computation  may  employ  one  million  particles, 
reach  steady  state  after  50,000  iterations,  and  continue  sampling  for  a  further  50,000  iterations.  On  a  modern 
desktop  computer,  such  a  simulation  should  take  about  3  hours. 
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While  the  ideas  behind  the  DSMC  technique  are  simple,  implementation  in  an  algorithm  takes  on  many 
different  forms.  Specific  DSMC  algorithms  have  been  developed  for  vector  computers  [10]  and  parallel 
computers  [11,12].  Bird  has  focused  work  on  customizing  the  algorithm  to  achieve  efficient  performance  on 
single  processor  machines  [13].  Some  of  the  most  commonly  used  DSMC  codes  are  summarized  in  the 
second  part  of  this  article  [8]  in  which  examples  of  their  application  to  various  atmospheric  entry  flows  are 
also  discussed. 

Having  provided  a  general  overview  of  the  basic  elements  of  the  DSMC  method,  in  the  following  sections 
some  of  the  physical  models  are  described  that  are  most  critical  to  the  application  of  the  DSMC  technique  to 
analyze  hypersonic  entry  flows. 


3.0  PHYSICAL  MODELS  OF  THE  DSMC  TECHNIQUE 

In  this  section,  the  most  commonly  employed  physical  models  are  reviewed  for  DSMC  computation  of 
hypersonic  entry  flows.  The  basic  ideas  are  described  for  simulation  of  a  number  of  physical  phenomena 
including  momentum  exchange,  internal  energy  relaxation,  chemical  reactions,  and  gas-surface  interactions. 
Where  possible,  examples  are  provided  of  attempts  to  validate  these  models  using  laboratory  data. 

3.1  Elastic  Momentum  Exchange 

In  the  absence  of  internal  energy  exchange  and  chemical  reactions,  an  elastic  collision  between  two  particles 
leads  only  to  changes  in  their  velocity  (or  momentum)  components.  The  frequency  of  such  interactions  is 
determined  by  the  collision  cross  section  and  a  number  of  such  models  have  been  developed  for  DSMC.  The 
most  widely  used  forms  are  the  Variable  Hard  Sphere  (VHS)  [14]  and  the  Variable  Soft  Sphere  (VSS)  [15]. 
For  hypersonic  flow,  the  VHS  model  is  considered  sufficiently  accurate,  for  which  the  cross  section  is  given 
as: 


cr  =  cr, 


f  \-2co 

g 


ref 


gref  J 


(3.1) 


where  (Jref  and  gref  are  reference  values,  and  (0  is  related  to  the  viscosity  temperature  exponent. 
Specifically,  it  is  assumed  that  the  gas  viscosity  is  described  by  a  simple  temperature  relation: 


V  =  /V 


\0.5+«) 


T 

\Aref  J 


(3.2) 


and  the  relationship  between  the  reference  parameters  is  provided  by  Bird  [4],  Values  of  these  reference 
parameters  for  many  common  species  are  provided  in  Bird  [4]  and  these  are  generally  obtained  by  comparison 
with  measured  or  computed  viscosity  data.  For  the  VHS  model,  isotropic  scattering  is  assumed  in  which  the 
unit  vector  of  the  post-collision  relative  velocity  is  assigned  at  random  on  the  unit  sphere.  The  VSS  model 
represents  an  improvement  over  VHS  in  that  it  allows  collision  parameters  to  be  determined  through 
comparison  with  both  viscosity  and  diffusivity  data.  The  VSS  cross  section  is  the  same  as  the  VHS  model, 
but  the  scattering  angle  is  given  by: 
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j  =  2cos  l< 


^Y°| 

\d) 


(3.3) 


where  a  is  determined  from  diffusivity  data,  b  is  the  distance  of  closest  approach  and  d  is  the  collision 
diameter.  Again,  values  of  a  for  common  gases  are  provided  by  Bird  [4],  Note  that  a=l  corresponds  to  the 
VHS  model. 


One  of  the  most  common  test  cases  for  evaluation  of  the  DSMC  technique  in  simulating  nonequilibrium  flows 
involving  only  elastic  momentum  exchange  are  normal  shock  waves  of  noble  gases.  Figure  la  shows  the 
density  profile  through  a  normal  shock  wave  at  Mach  9  for  argon.  Detailed  measurements  were  obtained 
using  an  electron  beam  technique  by  Alsmeyer  at  a  number  of  different  Mach  numbers  in  both  argon  and 
nitrogen  [16].  Included  in  Fig.  la  are  both  DSMC  and  CFD  results.  The  DSMC  computations  used  the  VHS 
model  while  the  CFD  results  solved  the  Navier-Stokes  equations  with  the  viscosity  given  by  the  same  VHS 
parameters  used  in  DSMC.  The  comparison  shows  that  the  DSMC  technique  is  able  to  reproduce  the 
measured  data  very  accurately  whereas  CFD  predicts  a  shock  wave  that  is  too  thin.  Figure  lb  shows  the 
reciprocal  shock  thickness  (a  measure  of  the  density  gradient  at  x  /  ’k  -  0)  for  all  the  argon  shock  waves 
investigated  by  Alsmeyer  [16].  Once  again,  it  is  clear  that  DSMC  provides  excellent  agreement  with  the 
measurements  for  all  conditions  considered.  CFD  consistently  predicts  an  inverse  shock  thickness  that  is  too 
large  (that  is,  shocks  that  are  too  thin)  for  all  Mach  numbers  above  about  1.50.  Another  compelling  validation 
of  the  capability  of  DSMC  in  simulating  nonequilibrium  hypersonic  flows  is  provided  by  the  comparisons 
shown  in  Figs.  2a  and  2b.  These  plots  are  taken  from  the  study  by  Pham-Van-Diep  et  al.  [17]  in  which 
velocity  distribution  functions  were  examined  inside  a  normal  shock  of  helium  at  Mach  25.  Figure  2a  shows 
the  distributions  in  the  front  of  the  shock  for  the  parallel  (circles  and  solid  line)  and  perpendicular  (triangles 
and  dashed  line)  velocity  components.  Symbols  represent  electron  beam  measurements  and  lines  represent 
DSMC  computations  that  employed  a  detailed  Maitland-Smith  collision  model.  The  horizontal  axis  in  these 
plots  is  reversed  so  that  the  high  freestream  parallel  velocity  component  appears  as  a  peak  towards  the  left  of 
the  figure.  The  parallel  velocity  distribution  shows  a  strong  nonequilibrium  profile  with  the  higher  velocity, 
lower  temperature  peak  on  the  left,  and  a  higher  temperature,  lower  velocity  peak  towards  the  middle.  The 
distribution  of  the  perpendicular  component  is  centered  on  zero  and  also  consists  of  two  distinct  populations 
from  the  low  temperature  freestream  and  the  high  temperature  post-shock  regions.  The  profiles  in  Fig.  2b  are 
obtained  further  downstream  towards  the  middle  of  the  shock  and  continue  to  show  strongly  nonequilibrium 
behavior.  Clearly,  the  Navier-Stokes  equations,  that  are  based  on  a  small  perturbation  from  equilibrium,  are 
not  able  to  accurately  model  such  phenomena.  The  excellent  agreement  with  measured  data  shown  in  these 
plots  is  one  of  the  strongest  illustrations  of  the  ability  of  the  DSMC  technique  to  reproduce  nonequilibrium 
flow  at  the  level  of  the  distribution  functions. 


The  main  limitation  of  the  VHS/VSS  collision  models  is  their  reliance  on  the  ability  to  model  viscosity  and 
diffusivity  as  a  simple  function  of  temperature.  Even  for  the  common  gases  for  which  VHS/VSS  parameters 
are  provided  by  Bird  [4],  the  viscosity  dependence  on  temperature  may  change  over  a  sufficiently  wide 
temperature  range.  For  illustration  of  the  types  of  problems  encountered  with  the  VHS/VSS  models,  Figs.  3a 
and  3b  show  collision  cross  sections  as  a  function  of  collision  energy  for  two  different  collisions  involving 
electrons.  For  such  interactions,  direct  measurements  of  cross  sections  are  available  in  the  literature.  Various 
attempts  to  fit  the  measured  data  using  the  VHS  model  are  shown  that  demonstrate  limitations  of  the  VHS 
model  for  describing  these  interactions.  Indeed,  a  negative  value  of  co  must  be  employed  in  VHS  in  order  to 
approach  the  measured  trends.  More  details  of  these  cross  section  evaluations  are  provided  in  [18]. 
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3.2  Rotational  Energy  Exchange 

The  DSMC  technique  usually  simulates  the  internal  energy  modes  of  atoms  and  molecules  by  assigning 
rotational,  vibrational,  and  electronic  energies  to  each  particle.  In  hypersonic  flows,  generally  the  electronic 
modes  are  ignored,  as  is  the  case  for  CFD  studies.  Analysis  of  jet  flows  including  electronic  energy  is 
described  for  example  in  [19].  We  focus  here  on  rotational  and  vibrational  energy  exchange. 


The  rotational  mode  is  usually  simulated  using  a  classical  physics  approach  in  which  the  rotational  energy,  £r, 
is  assumed  continuously  distributed  at  equilibrium  according  to  a  Boltzmann  distribution: 


f(£r)d£r  = 
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C  V72-1 

£ 


r(c/2) 


| \kTJ  x  ^  kT]  VkTj 
where  C,  is  the  number  of  rotational  degrees  of  freedom  (=2  for  a  diatomic  molecule;  =3  for  a  polyatomic 
molecule),  k  is  Boltzmann’s  constant,  and  T is  the  temperature. 
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(3.4) 


When  a  particle  representing  a  molecule  is  injected  into  a  DSMC  computation,  it  is  given  an  initial  rotational 
energy  sampled  from  Eq.  (3.4).  The  rotational  energy  of  the  particle  can  change  through  collisions  with  other 
particles  and  through  collisions  with  a  solid  surface  (see  section  3.7).  In  a  continuum  analysis  of  rotational 
energy  exchange,  the  rotational  relaxation  equation  is  usually  employed: 


der  _  er  -  er 

dt  Tr 


(3.5) 


* 

where  er  is  the  specific  rotational  energy,  er  is  the  equilibrium  value  at  temperature  T,  and  Tr  is  the 
rotational  relaxation  time.  The  equivalent  DSMC  procedure  involves  evaluating  a  probability  of  rotational 
energy  exchange  for  each  collision  followed  by  appropriate  energy  exchange  mechanics  for  those  collisions 
that  lead  to  rotational  relaxation.  The  average  probability  of  rotational  energy  exchange  is: 


1 

T,.V 


(3.6) 


where  Zrot  is  the  rotational  collision  number,  Tt  is  the  translational  relaxation  time  that  is  equal  to  the  inverse 
of  the  collision  frequency,  v.  Boyd  [20]  developed  the  following  instantaneous  rotational  energy  exchange 
probability  based  on  Parker’s  model  [21]  for  the  rotational  collision  number  and  the  VHS  collision  model: 


(3.7) 


V  J 

where  £tot  is  the  total  collision  energy  (the  sum  of  the  translational  collision  energy  and  the  rotational  energy), 
T  is  the  characteristic  temperature  of  the  intermolecular  potential,  and  (Zrot\  is  the  limiting  value.  After 


evaluation  of  the  rotational  energy  exchange  probability,  a  random  number  is  used  to  decide  whether  the 
collision  leads  to  energy  exchange.  For  those  collisions  involving  rotational  energy  exchange,  the  Borgnakke- 
Larsen  (BL)  model  [22]  is  employed  to  assign  new  post-collision  rotational  energies.  The  BL  model  assumes 
local  thermodynamic  equilibrium  to  sample  the  fraction  of  the  total  collision  energy  due  to  rotation,  £rot  / £tot , 
from  the  following  expression: 
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(3.8) 


Once  the  new  rotational  energy  is  assigned,  the  remaining  energy  is  the  new  translational  energy  and  hence 
determines  the  new  post-collision  relative  velocity.  The  regular  DSMC  collision  mechanics  is  then  performed 
to  calculate  the  velocities  of  the  colliding  particles. 


Lumpkin  et  al.  [23]  noted  that  an  additional  correction  must  be  applied  to  any  DSMC  rotational  energy 
exchange  probability  in  order  to  make  Borgnakke-Larsen  exchange  mechanics  consistent  with  the  continuum 
rotational  relaxation  equation,  Eq.  (3.5).  The  form  of  the  correction  is: 


p  =  p 

particle  continuun 


v  4  -  loo) 


that  is  usually  close  to  a  factor  of  two. 


(3.9) 


While  the  rotational  energy  is  usually  simulated  in  the  classical  limit,  a  quantum  mechanical  approach  has  also 
been  developed  by  Boyd  [24]. 

A  detailed  set  of  experimental  measurements  of  hypersonic  normal  shock  waves  in  nitrogen  was  collected  by 
Robben  and  Talbot  [25].  Again,  an  optical  diagnostic  technique  was  employed  to  measure  both  the  density 
and  the  rotational  energy  distribution  function  through  the  shock  wave  for  a  number  of  different  Mach 
numbers.  Figure  4a  compares  DSMC  simulations  [24]  with  the  measurements  of  the  density  and  rotational 
temperature  profiles  at  a  Mach  number  of  12.9.  Clearly,  very  good  agreement  between  simulation  and 
measurement  is  obtained.  Figure  4b  shows  rotational  energy  distribution  functions  measured  at  four  different 
locations  in  this  same  shock  wave,  moving  from  upstream  to  downstream  through  (a)  to  (d).  Once  again,  the 
excellent  agreement  obtained  between  DSMC  and  experiment  at  the  level  of  the  distribution  function  is  one  of 
the  strongest  illustrations  of  the  ability  of  the  technique  to  accurately  simulate  nonequilibrium  phenomena. 
The  main  limitation  of  the  Robben  and  Talbot  data  is  that  it  was  collected  for  a  flow  with  a  total  temperature 
of  just  300  K.  There  is  a  strong  requirement  for  additional,  detailed  measured  data  sets  for  true  hypersonic, 
entry  conditions. 


3.3  Vibrational  Energy  Exchange 

The  simulation  of  vibrational  relaxation  follows  a  similar  procedure  to  that  for  rotation.  The  average 
probability  of  vibrational  energy  exchange  is  typically  evaluated  using  the  vibrational  relaxation  time  used  in 
hypersonic  CFD  models  due  to  Millikan  and  White  [26]  with  the  Park  high  temperature  correction  [27]: 


'  vib 


7 MW  "I"  TPark 


(3.10) 


In  order  to  accurately  reproduce  this  vibrational  relaxation  time  in  a  DSMC  computation,  due  to  its  complex 
temperature  dependence,  it  is  necessary  to  evaluate  a  collision  averaged  vibrational  exchange  probability  [28]. 
Unlike  rotational  relaxation,  a  quantum  mechanical  approach  is  almost  always  employed  for  simulation  of 
vibrational  energy  relaxation  in  hypersonic  flows.  A  quantized  vibrational  energy  exchange  model 
corresponding  to  the  classical  Borgnakke-Larsen  approach  was  formulated  by  Bergemann  and  Boyd  [29].  It 
involves  first  determining  the  maximum  vibrational  quantum  level  available  based  on  the  total  collision 
energy: 
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(3.11) 


where  stot  is  the  total  collision  energy  (the  sum  of  the  translational  collision  energy  and  the  vibrational 

energy),  | _ |  means  truncation,  and  6y  is  the  characteristic  temperature  for  vibration  of  the  molecule.  Then, 

the  post-collision  vibrational  quantum  number,  v ,  is  sampled  from: 
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(3.12) 


The  Lumpkin  et  al.  [23]  correction  factor  must  also  be  applied  to  the  vibrational  exchange  probability. 


More  detailed  vibrational  relaxation  models  for  DSMC  have  also  been  developed  and  applied  to  hypersonic 
flows.  For  example,  in  [30],  the  Multiple  Quantum-Step  Transition  (MQST)  model  was  developed  for  use  in 
DSMC.  Based  on  the  Forced  Harmonic  Oscillator  (FHO)  model  of  Kemer  [31],  the  probability  of  a  multi¬ 
quantum  activation  is  given  by: 


v,v+Av 
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(3.13) 


sl=aexp{-g*  lg)  (3.14) 

where  constants  a  and  g*  are  determined  from  molecular  parameters.  The  corresponding  deactivation 
probability  is  developed  in  [30]  based  on  detailed  balance.  The  MQST  model  was  compared  in  [30]  to  the 
Borgnakke -Larsen  (B-L)  approach  for  flow  of  nitrogen  at  a  velocity  of  5.1  km/s  over  a  10  cm  radius  sphere  at 
80  km  altitude.  Figure  5  compares  the  computed  vibrational  energy  distribution  functions  at  the  point  of 
maximum  vibrational  temperature  in  the  shock  layer.  The  Boltzmann  distribution  corresponds  to  the 
maximum  vibrational  temperature.  Clearly,  both  the  phenomenological  B-L  model  and  the  more  detailed 
MQST  model  predict  a  significant  degree  of  vibrational  nonequilibrium.  It  is  interesting  to  note  that  the  B-L 
model  predicts  a  higher  level  of  nonequilibrium. 

There  are  no  measurements  in  the  literature  of  vibrational  energy  distribution  functions  in  hypersonic  flows 
that  can  be  used  to  validate  the  DSMC  vibrational  relaxation  models.  One  of  the  most  useful  sets  of  data  was 
obtained  by  Sharma  [32]  in  a  shock-tube  for  flow  of  N2  at  a  velocity  of  6.2  km/s.  Using  spectroscopy,  both 
the  rotational  and  vibrational  temperatures  of  N2  were  inferred  at  two  different  points  in  the  shock  wave  as 
shown  in  Fig.  6.  The  DSMC  computations  [28]  shown  in  Fig.  6  employed  the  B-L  approach  and  show 
reasonable  agreement  with  the  measured  data.  Further  more  detailed  experiments  of  this  kind  are  required  to 
more  completely  validate  the  DSMC  simulation  approach  for  vibrational  relaxation. 

3.4  Chemical  Reactions 

The  most  commonly  used  DSMC  chemistry  model  is  the  Total  Collision  Energy  (TCE)  model  of  Bird  [4], 
This  model  is  based  on  a  modified  Arrhenius  rate  coefficient  of  the  form: 
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C  =  aTh  exp 


kTJ 


(3.15) 


where  a  and  b  are  constants,  and  Eact  is  the  activation  energy  of  the  reaction.  By  integrating  over  the 
equilibrium  distribution  function  for  the  total  collision  energy,  it  may  be  shown  that  the  form  of  the  reaction 
probability  consistent  with  Eq.  (3.15)  for  the  VHS  collision  model  is  given  by: 


,  .  i 

(F  -F  )  ^  +2 

Ptce  =  ^,0‘  (3-16) 

where  stot  is  the  total  collision  energy  of  all  modes  of  both  particles  participating  in  the  collision,  and  constant 
A  depends  on  the  Arrhenius  parameters  and  molecular  constants.  The  TCE  model  was  extended  to  include  the 
important  physical  phenomenon  of  vibration-dissociation  coupling  by  Haas  and  Boyd  [33]  in  the  Vibrationally 
Favored  Dissociation  (VFD)  model.  The  VFD  model  makes  it  possible  to  increase  the  dissociation 
probability  of  particles  having  large  vibrational  energy: 
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Values  of  the  VFD  parameters  for  air  molecules  have  been  determined  through  comparison  with  experimental 
data.  For  example,  Fig.  7  shows  comparisons  [33]  for  dissociation  incubation  distances  between  DSMC 
predictions  and  data  measured  by  Hornung  [34]  in  hypersonic  flows  of  N2. 


Measurements  of  reaction  cross  sections  of  interest  in  hypersonic  air  flows  are  not  generally  available,  except 
for  some  reactions  involving  electrons.  As  ab  initio  computational  chemistry  techniques  mature,  there  is  the 
hope  in  the  future  that  detailed  computed  data  bases  can  be  used  to  help  develop  more  accurate  DSMC 
chemistry  models.  One  example  of  such  a  database  was  constructed  using  a  Quasi-Classical  Trajectory  (QCT) 
method  by  Bose  and  Candler  [35]  for  the  Zeldovich  exchange  reaction: 

N2  +  O^NO+N  (3.18) 


This  QCT  database  was  employed  [36]  to  perform  a  detailed  evaluation  of  the  TCE  model  for  this  particular 
reaction.  Figure  8a  shows  reaction  cross  sections  as  a  function  of  translational  collision  energy,  at  a  particular 
rotational  energy  level  (J=64)  of  the  reactant,  N2,  for  three  different  reactant  vibrational  levels  (v=0,  7,  13). 
Clearly,  the  TCE  model  is  not  at  all  accurate  for  this  reaction  in  terms  of  collision  cross  section.  This  poor 
comparison  motivated  the  development  of  another  DSMC  chemistry  model  that  allows  favoring  from  each  of 
the  translational,  rotational  and  vibrational  energy  modes.  Termed  the  Generalized  Collision  Energy  (GCE) 
model,  the  reaction  probability  is  given  by  [36]: 


(s  —  £  V 
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(3.19) 


Figure  8b  provides  comparisons  between  the  QCT  data  and  the  GCE  reaction  cross  sections.  While  the 
comparisons  are  far  from  perfect,  they  represent  a  significant  improvement  over  the  TCE  model.  For  this 
particular  reaction,  the  GCE  model  parameter  values  were:  a=0.2,  [3 =-0.5,  4>=0.3. 
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The  detailed  QCT  database  makes  it  possible  to  evaluate  several  aspects  of  DSMC  chemistry  modeling.  The 
second  aspect  concerns  the  energy  distribution  of  the  particles  selected  for  reaction.  Figure  9a  shows  the 
vibrational  energy  distribution  function  of  the  reacting  N2  molecules  under  a  thermal  equilibrium  condition 
where  all  mode  temperatures  are  at  5,000  K.  Clearly,  the  GCE  model  provides  almost  perfect  agreement  with 
the  QCT  data.  A  similar  level  of  agreement  is  shown  in  Fig.  9b  for  a  strongly  nonequilibrium  condition 
similar  to  that  expected  in  hypersonic  entry  flow  at  high  altitude. 

When  it  is  determined  in  a  DSMC  computation  that  a  chemical  reaction  occurs,  the  activation  energy  is 
removed  from  the  total  collision  energy,  then  Borgnakke-Larsen  procedures  are  used  to  distribute  the 
remaining  energy  among  the  energy  modes  of  the  product  particles.  Once  again,  the  QCT  database  allows  this 
aspect  of  DSMC  chemistry  modeling  to  be  assessed.  Figures  10a  and  10b  show  the  vibrational  energy 
distributions  of  the  product  NO  molecules  for  the  same  two  conditions  shown  in  Figs.  9a  and  9b.  Once  again, 
the  GCE  model  provides  remarkably  accurate  simulation  of  the  detailed  phenomena  captured  by  the  QCT 
analysis.  It  is  planned  in  the  near  future  to  conduct  further  comparisons  of  this  nature  for  important 
hypersonic  flow  reactions  such  as  dissociation  and  ionization. 

Recently,  details  of  simulating  backward  chemical  rate  processes  with  the  DSMC  technique  have  been 
discussed  by  Boyd  [37].  An  important  issue  here  is  that  the  TCE  model,  and  its  derivatives  like  VFD  and 
GCE,  all  depend  on  the  use  of  a  rate  coefficient  expressed  in  Arrhenius  form.  Backward  rate  coefficients  are 
generally  evaluated  as  the  quotient  of  the  forward  rate  coefficient  and  the  equilibrium  constant  that  is  usually  a 
highly  complex  function  of  temperature.  Thus,  it  is  not  generally  possible  to  express  the  backward  rate 
coefficient  in  the  simple  Arrhenius  form.  This  problem  is  either  addressed  by  performing  a  best-fit  of  the  true 
backward  rate  coefficient  to  an  Arrhenius  form  over  a  temperature  range  of  interest,  or  by  calculating  the 
temperature  in  each  cell  of  the  DSMC  computation  in  order  to  evaluate  the  equilibrium  constant  exactly  [37], 

While  the  DSMC  method  is  based  on  the  premise  of  a  dilute  gas  for  which  three  body  collisions  are  ignored,  it 
is  sometimes  important  to  include  recombination  in  DSMC  calculations  and  models  for  such  reactions  are 
presented  in  [4]  and  [28]. 

Several  other  DSMC  chemistry  models  relevant  to  hypersonic  entry  flows  have  been  developed  including  the 
maximum  entropy  model  [38],  the  weak  vibrational  bias  model  [39],  the  threshold  line  model  [40],  and  the 
model  of  Rebick  and  Levine  [41].  Many  of  the  models  are  reviewed  and  evaluated  in  [42]  and  [43]. 

3.5  Charged  Species 

Ions  and  electrons  are  formed  in  sufficiently  energetic  hypersonic  flows  first  through  associative  ionization 
and  then  through  electron  impact  ionization.  While  these  reactions  can  be  simulated  using  the  TCE  chemistry 
model  with  DSMC,  the  presence  of  electrons  in  a  DSMC  computation  presents  some  special  challenges. 
Specifically,  due  to  their  very  low  mass  relative  to  other  species:  (1)  the  thermal  velocity  of  an  electron  is 
orders  of  magnitude  higher  than  for  other  species;  and  so  (2)  the  collision  frequency  of  electrons  is  orders  of 
magnitude  higher  than  for  other  species.  If  an  electron  was  not  charged,  then  the  first  problem  would  simply 
require  use  of  a  much  smaller  time  step  At  than  would  be  required  otherwise.  However,  electrostatic  attraction 
means  that  electrons  and  ions  interact  with  one  another  such  that  electron  diffusion  is  reduced  and  ion 
diffusion  is  slightly  increased.  A  common  model  to  describe  this  behavior  is  ambi-polar  diffusion  in  which  it 
is  assumed  that  ions  and  electrons  diffuse  at  the  same  rate.  Bird  [44]  first  introduced  a  DSMC  model  for 
ambi-polar  diffusion  in  which  every  electron  particle  was  tied  directly  to  the  ion  it  was  born  with.  These  pairs 
of  charged  particles  then  move  throughout  the  flow  domain  based  on  the  velocity  components  of  the  ion 
particle.  This  method  is  reasonably  successful,  but  it  is  difficult  to  implement  and  has  poor  performance  at 
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high  ionization  levels.  Carlson  and  Hassan  [45]  introduced  a  scheme  in  which  electric  fields  are  evaluated 
based  on  averaged  charged  particle  properties.  The  charged  particles  are  moved  in  these  fields  with  the 
electrons  processed  at  a  significantly  smaller  time  step  than  the  ions.  An  alternate  approach  was  proposed  by 
Boyd  [46]  in  which  electron  particles  are  moved  throughout  the  domain  based  on  cell-averaged  ion  velocities. 
The  electron  and  ion  particles  are  no  longer  tied  together  explicitly  making  the  method  much  easier  to 
implement  and  the  approach  is  found  to  be  generally  more  robust  than  Bird’s  technique.  Using  either  of  the 
Bird  [44]  or  Boyd  [46]  approaches  means  that  electron  particles  are  moved  at  the  time  scale  of  the  heavier 
species  and  so  there  is  no  need  to  reduce  the  simulation  time  step. 

The  second  issue  faced  in  simulating  electrons,  related  to  their  significantly  higher  collision  frequency,  must 
also  be  addressed.  The  obvious  choices  are  to  either:  (1)  reduce  the  global  time  step;  (2)  allow  electron 
particles  to  collide  more  than  once  over  each  iteration;  or  (3)  perform  sub-cycling  of  collisions.  Sub-cycling 
involves  calling  the  collision  subroutine  several  times  over  each  movement  iteration  so  that  the  number  of 
collision  pairs  to  be  tested  is  evaluated  several  times  using  a  sub-cycling  time  step  Atc  that  is  smaller  than  the 
global  time  step  At. 

3.6  Radiation 

Radiation  is  of  interest  in  hypersonic  flows  in  terms  of  the  emission  signature  and  at  very  high  entry  speeds 
for  the  radiative  component  of  vehicle  heating.  Emission  signatures  are  usually  simulated  decoupled  from  the 
flow  field,  and  examples  are  discussed  of  such  analysis  in  the  second  part  of  this  article  [8].  There  have  been 
several  DSMC  studies  on  emission  signatures  where  the  excited  states  of  interest  were  simulated  directly  as 
additional  chemical  species  [47,48].  A  key  issue  here  is  the  availability,  or  lack  of  it,  of  accurate  excitation 
rate  coefficients.  For  estimation  of  radiative  heating.  Bird  first  included  thermal  radiation  effects  in  DSMC 
for  analysis  of  an  aero-assist  orbital  transfer  vehicle  [49,43]  and  his  ideas  were  extended  by  Carlson  and 
Hassan  [50].  These  models  essentially  represent  an  extension  of  the  rotational  and  vibrational  relaxation 
models  using  Borgnakke-Larsen  energy  disposal.  As  a  phenomenological  approach,  these  models  performed 
reasonably  well,  but  this  is  an  area  where  further  research  is  warranted. 

3.7  Gas-Surface  Interaction 

The  most  important  outcome  from  most  DSMC  analyses  of  hypersonic  entry  flows  is  the  determination  of  the 
properties  at  the  vehicle  surface  and  in  particular  the  aerodynamic  forces  and  moments,  and  the  convective 
heat  transfer.  The  surface  properties  are  of  course  very  sensitive  to  the  model  assumed  in  DSMC  for  gas- 
surface  interaction.  The  most  common  gas-surface  interaction  model  used  in  DSMC  is  fully  diffuse  reflection 
in  which  a  particle  reflects  from  the  surface  with  new  velocity  components  that  are  sampled  from  Maxwellian 
distributions  characterized  by  the  wall  temperature  (note  that  the  velocity  component  normal  to  the  surface  is 
sampled  from  a  biased  Maxwellian  distribution).  In  the  diffuse  model,  the  particle’s  internal  energies  are  also 
sampled  from  the  appropriate  equilibrium  distribution,  such  as  Eq.  (3.4)  for  rotation,  using  the  wall 
temperature.  The  opposite  limit  to  diffuse  reflection  is  specular  reflection  in  which  the  only  change  to  the 
particle’s  properties  is  its  velocity  component  normal  to  the  surface  that  is  simply  reversed  in  sign.  Many 
DSMC  computations  use  an  accommodation  coefficient,  a,  to  simulate  a  combination  of  diffuse  and  specular 
reflections  such  that  a=l  is  fully  diffuse,  and  a=0  is  fully  specular.  Real  engineering  surfaces  generally 
require  a  value  in  the  range  of  a=0.8-0.9. 

Figure  1 1  shows  a  comparison  of  measured  [51]  and  computed  distributions  of  argon  atoms  reflecting  from  a 
platinum  surface.  Clearly,  the  measured  pattern  is  not  reproduced  by  either  of  the  specular  or  diffuse  reflection 
models  (and  indeed  cannot  be  reproduced  by  any  combination  of  the  two  models).  This  type  of  comparison 
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led  to  the  development  of  more  sophisticated  gas-surface  interaction  models  for  DSMC,  for  example  the 
Cercignani-Lampis-Lord  (CLL)  model  [52],  Such  models  tend  to  have  a  stronger  theoretical  basis,  like  using 
a  reciprocity  relation,  and  offer  more  control  through  use  of  additional  parameters.  Figure  12a  shows  DSMC 
computed  particle  reflection  distributions  for  Mach  10  flow  of  N2  over  a  flat  plate  [53]  for  a  variety  of 
different  gas-surface  interaction  models.  Fully  diffuse  reflection  gives  the  cosine  distribution.  Fully  specular 
reflection  gives  a  distribution  that  is  almost  tangent  to  the  surface.  Through  variation  of  the  accommodation 
coefficients  aM  and  aT  (the  parameters  in  the  CLL  model),  a  wide  range  of  reflected  distributions  can  be 
generated  that  lie  between  the  ideal  limits  of  diffuse  and  specular  reflection.  Figure  12b  compares  profiles  of 
the  velocity  component  parallel  to  the  flat  plate  measured  by  Cecil  and  McDaniel  [54]  using  Planar  Laser 
Induced  Fluorescence  (PLIF)  and  a  number  of  DSMC  computations  using  diffuse  reflection  and  the  CLL 
model  with  a  range  of  model  parameters.  Note,  at  the  surface,  that  both  measurements  and  computations  show 
a  significant  level  of  velocity  slip.  Such  comparisons  allow  identification  of  appropriate  parameter  values  for 
use  in  the  CLL  model.  However,  the  use  of  such  models  is  relatively  limited  due  to  the  lack  of  this  type  of 
basic  information  to  identify  parameter  values  for  real  systems  of  interest.  This  is  perhaps  another  area  where 
computational  chemistry  can  provide  data  bases  that  can  be  used  to  build  more  advanced  DSMC  physical 
models. 

There  are  two  important  phenomena  arising  from  gas -surface  interaction  under  high  Knudsen  number 
conditions:  velocity  slip  and  temperature  jump.  The  relatively  low  number  of  collisions  experienced  by  a  gas 
at  high  Kn  means  that  the  average  velocity  at  the  wall  has  a  finite  value,  even  for  a  surface  with  fully  diffuse 
reflection.  This  phenomenon  reduces  shear  stress  and  may  affect  separation.  Similarly,  due  to  the  low 
collision  rate,  the  temperature  of  the  gas  at  the  wall  is  not  equilibrated  with  the  surface.  In  a  hypersonic  flow 
where  the  wall  temperature  is  cooler  than  the  gas,  this  phenomenon  leads  to  a  reduction  in  heat  transfer.  These 
effects  are  naturally  included  in  a  DSMC  simulation  using  the  diffuse  reflection  and  CLL  models  whereas  the 
usual  approach  for  CFD  is  to  assume  no  slip  and  no  temperature  jump  at  a  wall.  The  omission  of  these  high 
Knudsen  number  surface  phenomena  in  CFD  partially  explains  some  of  the  differences  noted  in  the  detailed 
comparisons  under  hypersonic  flow  conditions  of  DSMC  and  CFD  reported  in  [6,7] .  One  approach  to  try  and 
extend  the  usefulness  of  CFD  into  the  high  Knudsen  number  range  is  to  employ  velocity  slip  and  temperature 
jump  models,  see  for  example  [7].  However,  while  some  of  these  models  do  improve  the  agreement  between 
CFD  and  DSMC  results  for  surface  quantities,  it  is  not  always  achieved  with  a  corresponding  improvement  in 
the  comparisons  of  the  flow  properties.  This  situation  again  illustrates  the  need  to  perform  non-continuum 
computations  of  high  Knudsen  number  flows  using  kinetic  methods  such  as  the  DSMC  technique. 


4.0  SUMMARY 

The  direct  simulation  Monte  Carlo  method  (DSMC)  has  evolved  over  more  than  40  years  into  a  powerful 
analysis  tool  for  computation  of  kinetic  nonequilibrium  hypersonic  entry  flows.  The  heart  of  the  technique  is 
its  detailed  treatment  of  collisional  phenomena  including  momentum  exchange,  relaxation  of  internal  energy 
modes,  chemistry,  radiation,  and  gas-surface  interaction.  In  this  article,  the  current  status  of  DSMC  models 
for  simulating  these  physical  phenomena  has  been  reviewed.  It  is  demonstrated  that  the  DSMC  technique  is 
able  to  simulate  accurately  strong  nonequilibrium  phenomena  generated  inside  strong  shock  waves  of  noble 
gases.  Detailed  validation  of  rotational  energy  relaxation  models  is  also  demonstrated  at  the  level  of  rotational 
energy  distribution  functions  measured  inside  a  strong  shock  wave.  However,  detailed  measurements  of  such 
phenomena  at  the  high  temperature  conditions  of  hypersonic  entry  flows  are  still  missing.  Both 
phenomenological  and  detailed,  quantum  transition  DSMC  models  for  vibrational  relaxation  were  discussed. 
There  are  presently  no  measurements  of  vibrational  energy  distributions  in  hypersonic  flows  that  allow 
evaluation  of  such  models.  Several  different  DSMC  chemistry  models  are  described  and  two  of  them 


RTO-EN-AVT-162 


16(1)  - 13 


Direct  Simulation  Monte  Carlo  for  Atmospheric  Entry 
1.  Theoretical  Basis  and  Physical  Models 


e-t&rlA-'irtvi  Artt/rtkr 


evaluated  using  detailed  information  obtained  from  ab  initio  computational  chemistry  analysis.  In  the  absence 
of  detailed  measurements  of  reaction  cross  sections,  such  analyses  appear  very  promising  to  help  in  the 
construction  and  evaluation  of  detailed  DSMC  thermochemistry  models.  Descriptions  were  also  provided  of 
the  present  status  for  DSMC  computation  of  charged  species  (ions,  electrons),  and  radiation.  These  are  areas 
where  further  work  is  required.  Finally,  the  status  for  DSMC  computation  of  gas -surface  interaction  was 
reviewed.  It  is  shown  that  the  idealized  models  of  diffuse  and  specular  reflection  are  not  able  to  accurately 
reproduce  reflections  measured  under  hypersonic  conditions.  While  more  sophisticated  gas-surface  interaction 
models  have  been  developed,  their  use  is  limited  due  to  the  difficulty  in  determining  constants  in  the  models. 
Again,  this  appears  to  be  an  area  where  ab  initio  computational  analysis  may  be  useful  for  building  improved 
DSMC  modeling  capabilities. 
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Figure  la:  Profiles  of  normalized  density  through  a  Mach  9  normal  shock  wave  of  argon: 

measurements  from  [16]. 


Figure  1b:  Reciprocal  shock  thickness  for  normal  shock  waves  of  argon:  measurements  from  [16]. 
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Figure  2a:  Velocity  distribution  functions  in  the  front  a  Mach  25  normal  shock  of  helium:  symbols  = 

experiments;  lines  =  DSMC  [17]. 
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Figure  2b:  Velocity  distribution  functions  in  the  middle  of  a  Mach  25  normal  shock  of  helium: 
symbols  =  experiments;  lines  =  DSMC  [17]. 
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Figure  3a:  Electron-oxygen  atom  collision  cross  sections  [18]. 


Figure  3b:  Electron-nitrogen  molecule  collision  cross  sections  [18]. 
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Figure  4a:  Profiles  of  density  and  rotational  temperature  in  a  Mach  12.9  normal  shock  of  N2  [24], 


Figure  4b:  Rotational  energy  distribution  functions  at  four 
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Figure  5:  Vibrational  energy  distribution  functions  in  the  shock  layer  of  N2  flow  over  a  sphere  at  5.1 

km/s  and  80  km  altitude  [30]. 
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Figure  6:  Temperature  profiles  in  a  Mach  17.8  normal  shock  wave  of  N2  [28]. 
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Figure  7:  Dissociation  incubation  distance  as  a  function  of  atom  mass  fraction  [33]. 
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Figure  9a:  Vibrational  energy  distribution  function  of  reactants  in  nitric  oxide  formation:  Tt=Tr=Tv=5,000K  [36]. 
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Figure  9b:  Vibrational  energy  distribution  function  of  reactants  in  nitric  oxide  formation:  Tt=14,000K, 

Tr=5,000K,  Tv=1 ,000K  [36]. 


Figure  10a:  Vibrational  energy  distribution  function  of  products  in  nitric  oxide  formation:  Tt=Tr=Tv=5,000K  [36]. 
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Figure  10b:  Vibrational  energy  distribution  function  of  products  in  nitric  oxide  formation: 
Tt=1 4,000K,  Tr=5,000K,  Tv=1 ,000K  [36]. 


Figure  11 :  Distributions  of  Ar  scattering  from  platinum:  measurements  from  [51]. 
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Figure  12a:  DSMC  computed  scattering  distributions  from  a  DSMC  simulation  of  hypersonic  flow  of 

N2  over  a  flat  plate  [53]. 
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